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We propose an experiment to explore the magnetic phase transitions in interacting fermionic 
Hubbard systems, and describe how to obtain the ferromagnetic phase diagram of itinerant electron 
systems from these observations. In addition signatures of ferromagnetic correlations in the observed 
ground states are found: for large trap radii (trap radius Rt > 4, in units of coherence length £), 
ground states are topological in nature — a "skyrmion" in 2D, and a "hedgehog" in 3D. 
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I. INTRODUCTION 



The simplest, and best studied model of itinerant fer- 
romagnetism (FM) is the Hubbard modet^. Shortly after 
its introduction, Nagaoka and Thouless proved^ 3 , that 
for an interaction U of infinite strength, doping one hole 
into a background of spins leads to FM. Yet for more 
than forty years, the fate of the Thouless-Nagaoka phase 
in the Hubbard model at finite doping and finite inter- 
action has not been fully resolved. Existing studies in- 
clude various perturbative 4 -^, variational 7 -^"^, slave 
bosonii^, Quantum Monte-Carlo 1 ^!! and DMFT-i£ cal- 
culations. For example, there is no consensus as to the 
values of the critical doping S cr (« 0.19 — 0.49) below 
which FM occurs and critical interaction U cr (w 63 — 77.7 
in units of the hopping amplitude t). Moreover, the var- 
ious approaches do not agree on the nature of the tran- 
sition, first order claimed by som o 4 ' 5 ' 11 second order by 
othcrs 1 ^ 4 -^. 

Experiments on optical lattices, which allow a tunable 
control of model parameters, offer an interesting opportu- 
nity to address these long-standing open questions. This 
setting was used to study correlations in cold bosonic 
systems experimentall y 16 ' 1 7 , and could also be applied 
to superconductivity of fermions. The superconducting 
transition, however, requires a low entropy state cooled 
to temperatures far below degeneracy, which still poses a 
significant challenge. Ferromagnetism offers a particular 
advantage in this respect, since, at least for the Nagaoka 
phase at U — oo, the absence of any other scales sets 
the transition temperature to be a finite fraction of the 
bandwidth, and entropy at the transition is not a small 
parameter (of order unity per carrier) . A recent experi- 
ment on two-component fermionic 40 K in optical lattices 
has reported a Mott insulator with a maximum achieved 
U/t ~ 40, T/T F = 0.28 and initial densities of more 
than one electron per site 1 ^. A FM phase presumably 
resides near the observed Mott phase, and conservative 
estimates place it at slightly less than one particle per 
site, U/t ~ 100 and T/T F ~ O.lW 5 , The experimental 
parameters required in a setup such as 1 ^ are estimated: 
lattice depth Vq/E t ~ 13 — 19, depending on the scatter- 
ing length a s , the coherence length £ = n ~ 0.2/im, 



and the trap radius Rt ~ 1.6/xm, obtained from the laser 
waist radius of 160/xm. 

In this paper we explore the possibility of using cold 
fermion gases to study FM in both 2D and 3D systems. 
For large trap sizes (radius Rt > 4 coherence lengths 
£) with the constraint of zero total magnetization and 
filling factor less than unity everywhere in the trap, we 
find that the "skyrmion" configuration (Fig. [1] (a)) has 
the lowest energy in 2D, while in 3D, the lowest energy 
configuration is the "hedgehog" (Fig. [T] (b)). Observing 
these ground states provides a simple way to map out the 
FM phase diagram, which would be the first result of its 
kind. 



(a) 



(b) 





FIG. 1: (Color online) (a)Skyrmion and (b) Hedgehog config- 
uration of spins. 



II. APPROACH 

Our approach to the problem is in marked contrast 
with previous studie o 19 ' 20 i 21 ' 22 , some of which address 
only the 2D case: rather than using cold atoms as a 
means to verify microscopic calculations, we propose us- 
ing them as a direct probe of itinerant FM in both 2D 
and 3D traps. Moreover, in the analysis of 2 ^ total magne- 
tization is not conserved; whereas, in most of the recent 



experiments 



23.24.25 



total magnetization is, in fact, con- 



served, as a consequence of isolation from the environ- 
ment and an absence of coupling between the effective 
spin degree of freedom and the rest of the system. 

Ferromagnetism has been studied extensively in the 
context of multi-component Bosc-Einstcin condensates 
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(spinor BECs). primarily in two dimensions, starting 
wit h 26 i 27 i 28 i 29 i 30 i 31 . However, magnetic order in these 
systems is a secondary consequence of condensation, and 
develops wherever the condensate forms. Fermions near 
the point of degeneracy, on the other hand, become FM 
in the absence of any other broken symmetry. Trapping 
potentials give rise to spatial density variations, so differ- 
ent regions of the fermionic fluid at or near degeneracy 
may not reach the Stoner instability^ 2 - simultaneously, re- 
sulting in separation between the FM and paramagnetic 
(PM) phases^, something intrinsically absent in spinor 
BECs and bulk materials. This phase separation ensures 
that the resulting ferromagnetic droplets have definite 
size, and this, in turn, allows one to measure the phase 
diagram, as discussed below. 

We consider an optical lattice containing fermionic 
atoms cooled to temperatures close to degeneracy. The 
role of spin is played by two states in the atoms' hy- 
pcrfinc multiplet. Hopping between lattice sites is de- 
termined by the overlap of the atomic wavefunctions on 
those sites, while the on-site interaction term is given 
by the s-wave scattering of "electrons" with opposite 
"spin"—. In a trapped gas of equal spin populations, 
ground states are constrained by vanishing magnetiza- 
tion: M = J dr M(r) = 0. 

What kind of ferromagnetic quantum state satisfies 
this condition? One possibility is the state \ip) = \S,Q) 
with maximal total M 2 = S(S + 1), but M z = 0. The 
expectation values of M. x and Ai y in this state are zero, 
ostensibly satisfying the requirement that the total mag- 
netization (ip\A4\ip) = 0. However, we can rule out such 
a state by observing that each component of the total 
magnetization operator is conserved, [H, M a ] = 0, so 
that the quadrupole operator Q a p = M. a M.f3 — \M 2 5 a f 3 
is also conserved: [H, Q a p] = 0. The state \ip) has a non- 
vanishing quadrupole moment, where {Q xx ) = (Q yy ) = 
—2{Q)zz = S(S + l)/6, and therefore cannot develop as 
a spin-isotropic paramagnetic state is cooled through the 
ferromagnetic transition. 



III. LANDAU GINZBURG TREATMENT 

To go further, we need to consider states of non- 
uniform magnetization in which M. — 0. To this end, 
we make a long- wavelength expansion of the total energy 
as a functional of the local magnetization. Such a long- 
wavelength treatment of the problem does not imply that 
the underlying nature of the system is classical. Indeed, 
long-wavelength actions of this sort have been used to 
great success in the analysis of one-dimensional quantum 
antifcrromagnets, where long-range order is completely 
absent^. 

The Landau-Ginzburg free energy functional which de- 
scribes long-wavelength configurations of the magnetiza- 



tion order parameter takes the form 

F LG =/drf|VM| 2 + ^(|ili| 2 + ^) 2 (1) 

Odd-power terms are ruled out by the time-reversal sym- 
metry of the free energy. Coefficients p and (3 are assumed 
to be positive and constant for simplicity, and the entire 
effect of the trap potentials is in the position dependence 
of a(r). We define R c such that a(r) < for r < R c and 
a(r) > for r > R c - We assume that the density ev- 
erywhere in the system is less than one electron per site, 
and conclude with a discussion of the remaining cases. 

Clearly, whatever magnetic moment is accumulated by 
any one region in the trap must be completely canceled 
by the rest of the trap. The configuration which connects 
any two regions could either be a domain wall or some 
form of a twist (either with vorticity or without). It is 
the competition between these three scenarios that de- 
termines the ground state. Due to global rotation invari- 
ance, it suffices to consider a single representative state 
from each class. 




FIG. 2: (Color online) Showing (a) domain wall and (b) 
single-axis twist configuration of spins 

In 3D the first candidate is a hedgehog (Fig. [T] (b)). 
Qualitatively, it consists of a core of radius ro, where 
magnetization is suppressed, a region of thickness L in 
which it continuously increases from to some fraction 
a | < 1 of the uniform value Mo, and the outer region, 
extending to the edge of the trap at radius Rt = R c - In 
a particular realization, the magnetization vector at any 
point is in the radial direction. A competing configura- 
tion has a core of radius r maximally polarized in one 
direction. A domain wall of thickness L connects the core 
to the outer region polarized at a fraction < a < 1 of 
the maximum Mq in the opposite direction (Fig. [2] (a)). 
The remaining possibility is a pure twist. The twist plane 
is globally fixed (a rotation of 9q about the z-axis), and 
the magnetization turns about the x-axis through an an- 
gle air, 1 < a < 2. The twist occupies a shell of thickness 
L outside of a maximally-polarized core of radius ro , and 
the remaining outer region is polarized along the final 
direction (Fig. [2]) (b)). Parameters a, ro, and L for each 
configuration are determined by numerically minimizing 
Flg for a fixed value of R c subject to the constraint of 
zero net magnetization. 



In 2D in addition to the domain wall, the single-axis 
twist and the hedgehog, which are just planar slices 
of their 3D cousins, there is another possibility: the 
skyrmion. This configuration has a maximally polarized 
core of radius r$ at the center, a twist through the an- 
gle air, 1/2 < a < 1 in a ring of thickness L, and an 
outer region, polarized along the final direction (Fig. [1] 
(a)). Unlike the single axis twist, the twist axis of the 
skyrmion is a function of position. The constraint for the 
2D candidates must be implemented explicitly, and the 
ground state is obtained the same way as in 3D. 

Before solving the full problem, some estimates are in 
order. In 2D a domain wall is approximately the suppres- 
sion of the order parameter in an annulus of thickness L 
and inner radius tq. Omitting a dimensionful prcfactor 
common to all configurations, and in the case of the do- 
main wall, ignoring the stiffness contribution, we find the 
free energy (£ = \j pjo. is the coherence length) 

ftw~^ (2) 

The energy is minimized by shrinking L to L oc £, how- 
ever, the magnetization constraint forces r oc R C1 and so 
Fuw ~ Rc- In the same units the skyrmion free energy 
has contributions from vorticity and from twisting, i.e. 

Jta-H^ + W^ + i) (3) 

The total magnetic moment of the skyrmion is zero for 
some a ~ 0(1). Fsk is minimized when both L, rg oc R c . 
As a result, the skyrmion free energy is roughly indepen- 
dent of the trap size. Analysis of the single axis twist 
is similar, however, implementing the constraint requires 
more twisting in the same volume than in the skyrmion 
case, and therefore this configuration should always be 
higher in energy. The 2D hedgehog free energy is esti- 
mated to be i*H g ~ 2 In + i(y) 2 , sum of vorticity and 
core energy, and one can show that Fkg ~ In R c /£. Thus, 
for large trap sizes we expect a skyrmion to form in the 
trap. Similar treatment of the 3D configurations shows 
that the hedgehog should be the ground state. Numerical 
computations confirm these estimates, as shown in Fig. 
El We find that in small traps in both 2D and 3D do- 
main walls are preferred, while for large traps in 3D 
[R c /£, > 7) the hedgehog has the lowest energy, and in 
2D (i? c /£ > 4), the skyrmion is the ground state. 

In the above derivation we tacitly assumed that that 
the whole system becomes FM, i.e. R c = Rt- This is 
close to the truth in 2D, since the density of states is 
practically flat, and the Stoner criterion for magnetic in- 
stability dictates that the entire system becomes FM at 
once, independent of the spatial density variations. In 
3D, however, the density of states is more complicated, 
thus certain parts of the trap will cross into the broken 
symmetry regime earlier than others. In the language 
of Landau-Ginzburg, it means that a changes sign from 
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FIG. 3: (Color online) (a) Free energy vs. trap size in two 
dimensions, • hedgehog, A skyrmion, * domain wall, ■ single 
axis twist. Insert: Log e (F) vs. R c /£- (b) Free energy vs. 
size of the correlated region in three dimensions for a fixed 
trap size, Rt = 75£, • hedgehog, ♦ domain wall, ■ single axis 
twist. Insert: Log e (F) vs. R c /(,- 



positive outside of some radius R c , to negative inside. In- 
corporating this into our analysis requires somehow sup- 
pressing the magnetization outside R c in the candidate 
states. This can be accomplished by half of a domain 
wall, e.g. M(r) diminishing to within some shell of 
thickness L\ of the critical radius R c and remaining to 
the edge of the trap Rt- Although the qualitative results 
of the calculations remain the same, the numerical prob- 
lem itself is modified: we minimize Flg with respect to 
parameters a, tq, L, and L\ for a given radius of the FM 
region (R c ), and a given radius of the trap (Rt)- Results 
for 3D in Fig. Improperly reflect these considerations. 

In the case when the trap has regions with density 
higher than one electron per site, we expect further 
phase separation. In these regions the chemical potential 
may enter the Mott gap and we would expect an anti- 
ferromagnetic Mott phase to set in with filling locked to 
one electron per site. In regions of even higher density, 
the chemical potential might cross the gap and emerge, 
once again, in an itinerant band. Thus, we expect the 
phase profile to have wedding cake structure^!. For traps 
with very high density in the middle, we expect a hole 
PM at the core, followed in turn by shells of hole FM, an 
anti-ferromagnetic Mott insulator, and an electron FM, 
finally ending with the electron PM exterior. For inter- 
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mediate densities, the Mott insulator could form a natu- 
ral core for the 3D hedgehog configurations. 



IV. DISCUSSION 

Experiments with ultra-cold atomic gases generally ex- 
ploit absorption imaging to detect the state of matter in 
the trap. For multi-component gases in-situ imaging of 
each individual component is possibl e) 23 ' 24 ' 25 ' 35 . In the 
case of FM systems this technique allows resolving the 
full magnetization integrated along the camera axis. Fig. 
2] shows in-situ images of integrated model magnetization 
for the hedgehog and the skyrmion. The critical radius 
can be clearly seen (and measured) in the 3D configura- 
tion, and it designates the phase boundary between the 
PM outside and the FM inside. In addition to the phase 
boundary and the profile of the magnetization, one could 
look for signatures of correlation in the shot nois o 36 ' 37 . 
The auto-correlator (M z (r)M 2 (r')) can be calculated us- 
ing the time-of-flight images, since it is related to an in- 
tegral of the frequency space magnetization correlator. 
Tracking the associated correlation length, which will di- 
verge as the system approaches the transition, is, po- 
tentially, a more sensitive tool and should indicate the 
proximity of the FM phase before the image itself shows 
the phase separation 3 -^. 

With the newly found signatures of correlation we can 
map out the phase diagram of an interacting fermion gas. 
The density of particles at the critical radius R c can be 
obtained from the absorption images. Properly normal- 
ized, e.g. the Mott insulator is exactly at half-filling, 
this density defines the critical doping of the system. In- 
teraction strength U/t ~ a s exp(2y / Vo/E r ), is given by 
the s-wave scattering length, a s and the lattice depth Vb, 
with a s and Vb (i.e. U and t) controlled independently^ 3 -. 
Thus scanning a s and Vq and measuring the critical dop- 
ing for each, we determine the phase diagram. As dis- 
cussed above, at higher densities the magnetization im- 
ages should exhibit an even richer shell structure. How- 
ever, by particle-hole symmetry all of the FM regions 
should have the same basic features, and each phase 
boundary present in the image, will provide a point on 
the full phase diagram. The finer details of the observed 
configuration can give us even more details about the pa- 
rameters of the effective Landau-Ginzburg description; 
and in turn, our observation of the Landau-Ginzburg pa- 



rameters can be used to evaluate the efficacy of various 
microscopic calculational schemes. Similarly, we could 
explore the vicinity of the FM quantum critical point, 
and compare the observations to the Hertz-Millis theory. 
In any real experiment a number of practical issues will 
play a role. In particular, the details of symmetry break- 
ing must be the same in successive experimental runs and 
in the different planes of a quasi-2D trap configuration. 
One way to ensure this is to impose a small position de- 
pendent external field. A weak interplanar coupling in 
quasi-2D lattice might also stabilize the long range order 
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FIG. 4: (Color online) False color plot of magnetization inte- 
grated in-situ, symmetry broken the along z-axis. Hedgehog 
imaged along the x-axis (left) and skyrmion imaged along the 
z-axis (right). 

in the trap. Our model neglects the fact that in experi- 
ments to date, strong interactions mix the various bands 
of the optical lattic e - 38 ' 39 . Elimination of this effect may 
require a lower interaction strength, which in turn will 
require a still lower temperature^ 3 -. 

In this paper we demonstrated that ground states of 
interacting fermions in both 2D and 3D traps are config- 
urations with non-trivial topologies. We've determined 
detectable signatures of these states, and proposed an 
experiment to map the phase diagram using these sig- 
natures. Such experiments offer the prospect of directly 
probing itinerant FM, studying the Stoncr criterion in 
strongly correlated magnets, and providing much needed 
experimental input in the debate on the phase diagram 
of the Hubbard model. 
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